#**********************************************************************************
# ccamFigs.r
# This file contains the code to plot figures for ccam.
#
# Author            : Chris Grandin
# Development Date  : December 2011 - September 2012
#
#**********************************************************************************

saveFig <- function(filename){
  # Save the currently plotted figure to disk
  # in the 'figDir' folder with the name
  # filename.plotType
  if(saveon){
    filename <- paste(figDir,filename,".",plotType,sep="")
    savePlot(filename,type=plotType)
    cat(paste("Saved figure ",filename,"...\n",sep=""))
  }
}

fig.fishingMortality <- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op <- par(no.readonly=T)
    plot(A$yr,A$ft[1,],type="l",xlab="Year",ylab="Fishing Mortality (/yr)",ylim=c(0,1.5*max(A$ft[1,])), xlim=c(A$yr[1],A$yr[length(A$yr)]), las=1)
    saveFig("fig.fishingMortality")
    par(op)
  }else{
    cat(paste("Error plotting 'Fishing mortality' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.catchFit <- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op <- par(no.readonly=T)
    plot(A$yr,A$ct[1,],type="l",xlab="Year",ylab="Catch",ylim=c(0,1.5*max(A$obs_ct,A$ct)), xlim=c(A$yr[1],A$yr[length(A$yr)]), las=1)
    points(A$yr,A$obs_ct[1,])
    saveFig("fig.catchFit")
    par(op)
  }else{
    cat(paste("Error plotting 'Catch fit' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.comm.age.residuals1	<- function(scenario){
#Need to change the residual calculations to multivariate logistic here.
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
     if(AgeLikelihood!=3){
		  op <- par(no.readonly=T)
		  par(mfcol=c(8,5),mar=c(0,0,0,1),oma=c(5,5,1,1), las=1)
		  par(xaxt="n",yaxt="n")
		  iyr=Acom_res[,1]
		  for(i in 1:length(iyr)){
			    r <- Acom_res[i,3:(nage-sage_c+3)]
			    iclr <- r
			    iclr[r>=0] <- 1
			    iclr[iclr!=1] <- 2
			    plot(sage_c:nage,r,type="h",ylim=c(-6,6),main="",col=iclr, las=1)
			    points(sage_c:nage,r,pch=19,col=iclr)
			    title(sub=iyr[i],line=-5)
		   }
			mtext("Residual",2,outer=T,las=3,line=2)
			mtext("Age",1,outer=T,line=1)
			saveFig("fig.comm.age.residuals1")
			par(op)
    }else{
      cat(paste("Error plotting 'Age Residuals' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Age Residuals' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.comm.age.residuals2 <-	function(scenario) {
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    if(AgeLikelihood!=3){
		  op <- par(no.readonly=T)
		  r <- Acom_res
			bubble.plot(r[,1],sage_c:nage,r[,3:(nage-sage_c+3)],scale=0.3,xlab="Year",ylab="Age",add=F,log.scale=T, las=1)
			par(mfcol=c(1,1))
			abline(a=2-1972,b=1,lty=2,col=3)
			abline(a=2-1975,b=1,lty=2,col=3)
			abline(a=2-1979,b=1,lty=2,col=3)
			abline(a=2-1982,b=1,lty=2,col=3,lwd=2)
			abline(a=2-1986,b=1,lty=2,col=3,lwd=2)
			abline(a=2-1989,b=1,lty=2,col=3)
			abline(a=2-1992,b=1,lty=2,col=3)
			abline(a=2-2001,b=1,lty=2,col=3,lwd=2)
		  saveFig("fig.comm.age.residuals2")
		  par(op)
    }else{
      cat(paste("Error plotting 'Age Residuals' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Age Residuals' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.comm.age.props	<-	function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    if(AgeLikelihood!=3){
      op <- par(no.readonly=T)
      par(mfcol=c(1,1))
      r <- Acom_obs

      bubble.plot(r[,1],sage_c:nage,r[,3:(nage-sage_c+3)],scale=0.3,xlab="Year",ylab="Age",add=F,log.scale=T, las=1)
      abline(a=2-1972,b=1,lty=2,col=3)
      abline(a=2-1975,b=1,lty=2,col=3)
      abline(a=2-1979,b=1,lty=2,col=3)
      abline(a=2-1982,b=1,lty=2,col=3,lwd=2)
      abline(a=2-1986,b=1,lty=2,col=3,lwd=2)
      abline(a=2-1989,b=1,lty=2,col=3)
      abline(a=2-1992,b=1,lty=2,col=3)
      abline(a=2-2001,b=1,lty=2,col=3,lwd=2)

	  	saveFig("fig.comm.age.props")
      par(op)
    }else{
      cat(paste("Error plotting 'Age Proportions' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Age Proportions' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.comm.age.props.fit	<-	function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    if(AgeLikelihood!=3){	
		  #Plots the collapsed conditional age-length keys to proportions at age in the survey
      op <- par(no.readonly=T)
      par(mfcol=c(8,5),mar=c(0,0,0,1),oma=c(5,5,1,1), las=1)
      par(xaxt="n",yaxt="n")

      iyr=Acom_obs[,1]
      for(i in 1:length(iyr)){
				mp <- barplot(Acom_obs[i,3:(nage-sage_c+3)],names.arg=paste(sage_c:nage),ylim=c(0.,0.8),main="", las=1)
				lines(mp,Acom_est[i,3:(nage-sage_c+3)],type="o",pch=20,cex=2,col=2)
				title(sub=iyr[i],line=-5)
      }
      mtext("Proportion",2,outer=T,las=3,line=2)
      mtext("Age",1,outer=T,line=1)

      saveFig("fig.comm.age.props.fit")
      par(op)
    }else{
      cat(paste("Error plotting 'Age Commercial Age Proportions fit' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Age Commercial Age Proportions fit' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.survey.age.residuals1	<-	function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    if(AgeLikelihood!=3){	
		  #Need to change the residual calculations to multivariate logistic here.
      op <- par(no.readonly=T)
		  par(mfcol=c(5,3),mar=c(3,2,2,1),oma=c(3,5,1,2), las=1)

		  iyr=Asurv_res[,1]
		  for(i in 1:length(iyr)){
        r <- Asurv_res[i,3:(nage-sage_s+3)]
        iclr <- r
        iclr[r>=0] <- 1
        iclr[iclr!=1] <- 2

        plot(sage_s:nage,r,type="h",ylim=c(-6,6),main=iyr[i],col=iclr, las=1)
        points(sage_s:nage,r,pch=19,col=iclr)
			}
			mtext("Residual",2,outer=T,las=3,line=2)
			mtext("Age",1,outer=T,line=1)

      saveFig("fig.survey.age.residuals1")
			par(op)
    }else{
      cat(paste("Error plotting 'Survey Age Residuals' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Survey Age Residuals' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.survey.age.residuals2 <-	function(scenario) {
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    if(AgeLikelihood!=3){
      op <- par(no.readonly=T)
      r <- Asurv_res
      bubble.plot(r[,1],sage_s:nage,r[,3:(nage-sage_s+3)],scale=0.3,xlab="Year",ylab="Age",add=F,log.scale=T, las=1)
      par(mfcol=c(1,1))
      abline(a=2-1986,b=1,lty=2,col=3,lwd=2)
      abline(a=2-1989,b=1,lty=2,col=3)
      abline(a=2-1992,b=1,lty=2,col=3)
      abline(a=2-2001,b=1,lty=2,col=3,lwd=2)
      saveFig("fig.survey.age.residuals2")
      par(op)
    }else{
      cat(paste("Error plotting 'Survey Age Residuals' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Survey Age Residuals' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.survey.age.props	<-	function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    if(AgeLikelihood!=3){
      op <- par(no.readonly=T)
      par(mfcol=c(1,1))
      r <- Asurv_obs
      bubble.plot(r[,1],sage_s:nage,r[,3:(nage-sage_s+3)],scale=0.3,xlab="Year",ylab="Age",add=F,log.scale=T, las=1)
      abline(a=2-1986,b=1,lty=2,col=3,lwd=2)
      abline(a=2-1989,b=1,lty=2,col=3)
      abline(a=2-1992,b=1,lty=2,col=3)
      abline(a=2-2001,b=1,lty=2,col=3,lwd=2)

      saveFig("fig.survey.age.props")
      par(op)
    }else{
      cat(paste("Error plotting 'Survey Age Proportions' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Survey Age Proportions' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.survey.age.props.fit	<-	function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    if(AgeLikelihood!=3){
		  #Plots the collapsed conditional age-length keys to proportions at age in the survey
      op <- par(no.readonly=T)
      par(mfcol=c(5,3),mar=c(3,2,2,1),oma=c(3,5,1,2), las=1)

      iyr=Asurv_obs[,1]
      for(i in 1:length(iyr)){
				mp <- barplot(Asurv_obs[i,3:(nage-sage_s+3)],names.arg=paste(sage_s:nage),ylim=c(0.,0.8),main=iyr[i], las=1)
				lines(mp,Asurv_est[i,3:(nage-sage_s+3)],type="o",pch=20,cex=2,col=2)
      }
      mtext("Proportion",2,outer=T,las=3,line=2)
      mtext("Age",1,outer=T,line=1)

      saveFig("fig.survey.age.props.fit")
      par(op)
    }else{
      cat(paste("Error plotting 'Survey Age Proportions fit' figure.  The model is ageless for scenario ",scenario,"\n\n"))
    }
  }else{
    cat(paste("Error plotting 'Survey Age Proportions fit' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

#RF NOTE - ADD AN MCMC VERSION OF THIS
fig.surveybiomass.fit	<-	function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op <- par(no.readonly=T)
    if(A$nits==1){
      plot(A$iyr,A$pit,type="l",xlab="Year",ylab="Survey biomass index",ylim=c(0,1.5*max(A$pit)), xlim=c(A$iyr[1],A$iyr[length(A$iyr)]), las=1)
      points(A$iyr,A$it)
      saveFig("fig.surveybiomass.fit")
    }
    if(A$nits>1){
      plot(A$iyr[1,], A$pit[1,], type="l",xlab="Year",ylab="Survey biomass index",ylim=c(0,1.5*max(A$pit[1,])), xlim=c(A$iyr[1,1],A$iyr[1,length(A$iyr[1,])]), las=1)
      points(A$iyr[1,], A$it[1,])
      windows()
      A$iyr2<-A$iyr[2,which(A$iyr[2,]>0)] #get rid of NAs
      A$pit2<-A$pit[2,which(A$pit[2,]>0)] #get rid of NAs
      A$it2<-A$it[2,which(A$it[2,]>0)] #get rid of NAs

      plot(A$iyr2,A$pit2,type="l",xlab="Year",ylab="Survey biomass index",ylim=c(0,1.5*max(A$pit2)), xlim=c(A$iyr2[1],A$iyr2[length(A$iyr2)]), las=1)
      points(A$iyr2,A$it2)
  		saveFig("fig.surveybiomass2.fit")
    }
    if(A$nits>2){
      windows()
      A$iyr3<-A$iyr[3,which(A$iyr[3,]>0)] #get rid of NAs
      A$pit3<-A$pit[3,which(A$pit[3,]>0)] #get rid of NAs
      A$it3<-A$it[3,which(A$it[3,]>0)] #get rid of NAs
      plot(A$iyr3,A$pit3,type="l",xlab="Year",ylab="Survey biomass index",ylim=c(0,1.5*max(A$pit3)), xlim=c(A$iyr3[1],A$iyr3[length(A$iyr3)]), las=1)
      points(A$iyr3,A$it3)
  		saveFig("fig.surveybiomass3.fit")
    }
    if(A$nits==4){
      windows()
      A$iyr4<-A$iyr[4,which(A$iyr[4,]>0)] #get rid of NAs
      A$pit4<-A$pit[4,which(A$pit[4,]>0)] #get rid of NAs
      A$it4<-A$it[4,which(A$it[4,]>0)] #get rid of NAs
      plot(A$iyr4,A$pit4,type="l",xlab="Year",ylab="Survey biomass index",ylim=c(0,1.5*max(A$pit4)), xlim=c(A$iyr4[1],A$iyr4[length(A$iyr4)]), las=1)
      points(A$iyr4,A$it4)
  		saveFig("fig.surveybiomass4.fit")
    }

    par(op)
  }else{
    cat(paste("Error plotting 'Survey Biomass fit' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }

}

fig.selectivity	<-	function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op <- par(no.readonly=T)
    sel<-exp(A$log_sel)
    sel[,1]<-log(sel[,1])
    vva<-sel[which(sel[,1]==1),]
    va<-vva[length(vva[,1]),2:length(vva[1,])]
    leg <- "Gear 1" #legend text
    cha <- 1 #plot character for legend
    plot(age, va,type="o",xlab="Age",ylab="Selectivity in final year",ylim=c(0,1), las=1)
    for(i in 2:ngear) {
      vva2<-sel[which(sel[,1]==i),]
      va2<-vva2[length(vva2[,1]),2:length(vva2[1,])]
      lines(age,va2,type="o",pch=17+i)
      leg <- c(leg,paste("Gear",i))
      cha <- c(cha, 17+i)
    }
    legend("bottomright",leg,lty=1,pch=cha,bty="n", cex=1.9)
    saveFig("fig.selectivities")
    par(op)
  }else{
    cat(paste("Error plotting 'Selectivity' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.phase <- function(yUpperLimit=2,scenario){
	#The phase plot of Bt/Bmsy vs (1-SPR)/(1-SPRmsy)
  # yUpperLimit is the upper limit of the y-axis
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    spr <- A$sprmsy_status

    sbstatus <- A$sbtmsy_status
    sbstatus <- sbstatus[1:length(A$yr)]

    maxX <- max(sbstatus)
    maxY <- yUpperLimit

    plot(sbstatus, spr, type="n", col=1, ylim=c(0, maxY), xlim=c(0, maxX), xlab="B/BTarget", ylab="(1-SPR)/(1-SPRTarget)")
	  #dum <- fried.egg(sbstatus, spr)

    lines(sbstatus, spr, type="o",col="blue")
    colVector <- vector(length=length(A$yr))
    colVector[1] <- "green"
    colVector[length(colVector)] <- "red"
    for(i in 2:(length(colVector)-1)){
      colVector[i] <- 0
    }
    points(sbstatus, spr, pch=20, col=colVector)
    abline(h=1, v=1, lty=3,col="red")

    saveFig("fig.phase")
    par(op)
  }else{
    cat(paste("Error plotting 'Selectivity' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.estimated.params.pairs <- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op <- par(no.readonly=T)
    ngr<-A$ngear #number of gears
    nits<-A$nits #number of survey gears
    npar <- A$npar #number of columns to plot starts out as number of estimated parameters

  	#get only burned in samples
    mc<-A$mc[Burn:nrow(A$mc),]

	  #don't want all columns -  want estimated pars (theta), ngear x ahat, and nits x q
	  #get column indices
    grcols<-vector(length=ngr)
    k=1
    for(i in 1:ngr){
      grcols[i]<-npar+k
      k=k+2 #skip ghat
    }

    nitcols<-vector(length=nits)
    for(i in 1:nits){
      nitcols[i]<-npar+k
      k=k+1
    }
    paircols<-c(1:npar,grcols,nitcols)

    j<-mc[,paircols]
    pairs(j, pch=".",upper.panel=panel.smooth,diag.panel=panel.hist, lower.panel=panel.smooth)
    saveFig("fig.estimated.params.pairs")
    par(op)
  }else{
    cat(paste("Error plotting 'Estimated Parameter Pairs' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}


fig.mcmc.priors.vs.posts <- function(exFactor=1.0,showEntirePrior=T,scenario){
	  # exFactor is a multiplier for the minimum and maximum xlims.
	  # qPriorFunction: 4 for gamma
	  # ghat=1 for survey selectivity, 2 for commercial selectivity
	  # showEntirePrior, if T then plot the entire prior function to its limits
	  #  and ignore posterior distribution limits
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
	  op	<- par(no.readonly=T)
	  ngr<-A$ngear #number of gears
	  nits<-A$nits #number of survey gears
	  npar <- A$npar #number of columns to plot starts out as number of estimated parameters
    nPar<-npar+nits+2*ngr
	  par(mfrow=c(4,4),omi=c(.1,.2,.2,.1), mai=c(.6,.5,.3,.1))  # show all parameters both fixed and estimated

	  # The values in the REPORT file for each of priorN are:
	  # 1. ival  = initial value
	  # 2. lb    = lower bound
	  # 3. ub    = upper bound
	  # 4. phz   = phase
	  # 5. prior = prior distribution funnction
	  #             0 = Uniform
	  #             1 = normal    (p1=mu,p2=sig)
	  #             2 = lognormal (p1=log(mu),p2=sig)
	  #             3 = beta      (p1=alpha,p2=beta)
	  #             4 = gamma     (p1=alpha,p2=beta)
	  # 6. p1 (defined by 5 above)
	  # 7. p2 (defined by 5 above)
	  functionNames <- c(dunif,dnorm,dlnorm,dbeta,dgamma)

  	###########################################################################
	  #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
	  #Get matrix of posteriors
	  #get only burned in samples
    mc<-A$mc[Burn:nrow(A$mc),]

	  #don't want all columns -  want estimated pars (theta), and nits x q
	  #get column indices
    grcols<-vector(length=2*ngr)
    k=1
    for(i in 1:(2*ngr)){
      grcols[i]<-npar+k
      k=k+1
    }

    nitcols<-vector(length=nits)
    for(i in 1:nits){
      nitcols[i]<-npar+k
      k=k+1
    }
    paircols<-c(1:npar,grcols,nitcols)
    j<-mc[,paircols]
    nm <- colnames(j)
    print(j[1:10,])

  	#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  	# Get MPD values 
	  mle=vector(length=nPar)
	  for(i in 1:npar){
      mle[i] <- A$theta[i]
	  }
    k=1
    for(i in 1:ngr){
	  	mle[npar+k] <- A$ahat[i]
	  	mle[npar+1+k] <- A$ghat[i]
	  	k=k+2
	  }
	  for(i in 1:nits){
      mle[npar+k] <- A$q[i]
	  }
	  print(mle)

	  #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
	  #Get priors
	  # the following are +1 because they are 0-based in the CTL file,
	  # and R's vectors are 1-based
	  fn<- vector(length=nPar)
	  mu<- vector(length=nPar)
	  sig<- vector(length=nPar)
	  for(i in 1:npar){
      fn[i]<-functionNames[A$theta_control[i,5]+1]
      mu[i]<-A$theta_control[i,6]
      sig[i]<-A$theta_control[i,7]
	  }
    k=1
	  for(i in 1:ngr){
      fn[npar+k]<-functionNames[1] #ahat (always uniform)
      if(A$ghat_pswitch[i]==1) fn[npar+1+k]<-functionNames[5] #gamma prior on ghat[i]
      if(A$ghat_pswitch[i]==0) fn[npar+1+k]<-functionNames[1] #uniform prior on ghat[i]

      mu[npar+k]<-0
      if(A$ghat_pswitch[i]==1) mu[npar+1+k]<-A$ghat_p1
      if(A$ghat_pswitch[i]==0) mu[npar+1+k]<-0.00000001

      sig[npar+k]<-A$age[length(A$age)]
      if(A$ghat_pswitch[i]==1) sig[npar+1+k]<-A$ghat_p2
      if(A$ghat_pswitch[i]==0) sig[npar+1+k]<-100

      k=k+2

    }
    for(i in 1:nits){
      fn[npar+k]<-functionNames[A$q_prior[i]+1]
      mu[npar+k]<-A$q_mu[i]
      sig[npar+k]<-A$q_sd[i]
	  }

	  for(param in 1:nPar){
			x <- list(p=j[,param],mu=mu[param],sig=sig[param],fn=fn[[param]],nm=nm[param],mle=mle[param])
			plot.marg(x,breaks=17,col="wheat",exFactor=exFactor,showEntirePrior)
	  }
	  saveFig("fig.mcmc.priors.vs.posts")
		par(op)
  }else{
    cat(paste("Error plotting 'Priors / Posteriors comparison' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

plot.marg <- function(xx,breaks="sturges",exFactor=1.0,showEntirePrior=T,...){
  #xx is a list(p=samples, mu=prior mean, s=prior varian, fn=prior distribution)
  # exFactor is a multiplier for the minimum and maximum xlims.
  # showEntirePrior, if T then plot the entire prior function to its limits
  #  and ignore posterior distribution limits
	ssNoPlot <- hist(xx$p,breaks=breaks,plot=F)
  xl <- seq(min(ssNoPlot$breaks),max(ssNoPlot$breaks),length=250)
  pd <- xx$fn(xl,xx$mu,xx$sig)
  z <- cbind(xl,pd)
  if(showEntirePrior){
    xlim <- c(min(xl),max(xl))
  }else{
    xlim <- c(min(ssNoPlot$breaks)/exFactor,max(ssNoPlot$breaks)*exFactor)
  }
  ss <- hist(xx$p,prob=T,breaks=breaks,main="",xlab=xx$nm,xlim=xlim,...)
  lines(xl,pd,col="green",lwd=2)     
  abline(v=xx$mle, lwd=2, lty=2, col=2)
}

fig.mcmc.trace <- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T) 
    par(mfrow=c(4, 4), las=1)
    mcmcData <- window(mcmc(A$mc[,1:16]), start=Burn, thin=Thin)
    for(param in 1:ncol(mcmcData)){
      par(mar=c(2,3,2,2))
      mcmcTrace <- as.matrix(mcmcData[,param])
      plot(mcmcTrace,main=colnames(mcmcData)[param],type="l",ylab="",xlab="",axes=F)
      box()
      at <- seq(0,end(mcmcData)-start(mcmcData),200)
      labels <- seq(start(mcmcData),end(mcmcData),200)
      axis(1,at=at,labels=labels)
      axis(2)
    }
    saveFig("fig.mcmc.trace")
    par(op)
  }else{
    cat(paste("Error plotting 'Traces' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.mcmc.density <- function(color=1,opacity="20",scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T) 
    par(mfrow=c(4, 4), las=1)
    mcmcData <- window(mcmc(A$mc[,1:16]), start=Burn, thin=Thin)
    for(param in 1:ncol(mcmcData)){
      par(mar=c(2,3,2,2))
      dens <- density(mcmcData[,param])
      plot(dens,main=colnames(mcmcData)[param])
      xx <- c(dens$x,rev(dens$x))
      yy <- c(rep(min(dens$y),length(dens$y)),rev(dens$y))
      shade <- getShade(color,opacity)
      polygon(xx,yy,density=NA,col=shade)
    }
    saveFig("fig.mcmc.density")
    par(op)
  }else{
    cat(paste("Error plotting 'Estimated Parameter Density' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.mcmc.autocor <- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    n <- 0
    par(mfrow=c(4, 4), las=1)
    for(i in 1:16){
      n <- n+1
      x <- window(mcmc(A$mc[,i]),start=Burn,thin=Thin)
      ac <- autocorr(x, lags = c(0, 1, 5, 10, 15,20,30,40,50), relative=TRUE)
      par(mar=c(2,3,2,2))
      autocorr.plot(x, lag.max=100,  auto.layout=F)
    }
    saveFig("fig.mcmc.autocor")
    par(op)
  }else{
    cat(paste("Error plotting 'Estimated Parameter Autocorrelation' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.mcmc.gelman <- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    x <- window(mcmc(A$mc[,1:12]),start=Burn+1,thin=Thin)
    gelman.plot(x,bin.width = 10,max.bins = 50,confidence = 0.95, transform = FALSE, autoburnin=TRUE, auto.layout = TRUE)
    saveFig("fig.mcmc.gelman")
    par(op)
  }else{
    cat(paste("Error plotting 'Estimated Parameter Gelman' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.mcmc.geweke <- function(frac1=0.1,frac2=0.5,nbins=20,pvalue=0.05,silent=F, useShades=F, scenario, ...){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op	<- par(no.readonly=T)
    on.exit(par(op))
    x <- window(mcmc(A$mc[,1:12]),start=Burn+1,thin=Thin)
    #x <- mcmc(read.table("gewekeExample.csv",sep=",",header=T)[,1:10])
    par(mfrow=c(4,3), las=1)
    ystart <- seq(from = start(x), to = (start(x) + end(x))/2, length = nbins)
    gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 
                 dimnames = c(ystart, varnames(x), chanames(x)))
    for (i in 1:length(ystart)) {
      geweke.out <- try(geweke.diag(window(x, start = ystart[i]),frac1 = frac1, frac2 = frac2), silent=silent)
      for (k in 1:nchain(x)){
        gcd[i, , k] <- geweke.out[[k]]
      }
    }

    climit <- qnorm(1 - pvalue/2)
    for (k in 1:nchain(x)){
      for (j in 1:nvar(x)) {
        par(mar=c(2,3,2,2))
        ylimit <- max(c(climit, abs(gcd[, j, k])))
        plot(ystart, gcd[, j, k], type = "p", xlab = "Iteration", 
             ylab = "Z-score", pch = 1, ylim = c(-ylimit, ylimit), col="blue", 
             ...)
        lines(ystart,gcd[,,1][,j],lwd=2, col="blue")
        if(useShades){
          xx <- c(ystart,rev(ystart))
          yy <- c(rep(100*max(gcd[,,1][,j]),length(gcd[,,1][,j])),rev(gcd[,,1][,j]))
          shade <- getShade("green","15")
          polygon(xx,yy,density=NA,col=shade)

          xx <- c(ystart,rev(ystart))
          yy <- c(rep(100*min(gcd[,,1][,j]),length(gcd[,,1][,j])),rev(gcd[,,1][,j]))
          shade <- getShade("blue","15")
          polygon(xx,yy,density=NA,col=shade)

        }
        abline(h = c(climit, -climit), lty = 2)
        if (nchain(x) > 1) {
          title(main = paste(varnames(x, allow.null = FALSE)[j], 
                  " (", chanames(x, allow.null = FALSE)[k], ")", 
                  sep = ""))
        }else{
          title(main = paste(varnames(x, allow.null = FALSE)[j], 
                  sep = ""))
        }
      }
    }
    saveFig("fig.mcmc.geweke")
    par(op)
  }else{
    cat(paste("Error plotting 'Estimated Parameter Geweke' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.variance.partitions <- function(scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[7]]){ # if MCMC results are loaded
    op <- par(no.readonly=T)
    rho <- A$mc$rho[Burn:nrow(A$mc)]
    varphi <- 1/A$mc$varphi[Burn:nrow(A$mc)]
    sig <- rho*varphi
    tau <- (1-rho)*varphi
    d <- cbind(vp=varphi,sig,tau)
    pairs(d, pch=".", upper.panel=NULL, gap=0)
    saveFig("fig.variance.partitions")
    par(op)
  }else{
    cat(paste("Error plotting 'Estimated Parameter Variance Partitions' figure.  There are no MCMC outputs loaded for scenario",scenario,"\n\n"))
  }
}

fig.time.varying.selectivity	<- function(gear=1,scenario){
  try(dev.off(),silent=T)
  if(opList[[scenario]][[6]]){ # if MPD results are loaded
    op	<- par(no.readonly=T)
    with(A, {
      plot.sel<-function(x, y, z, ...){
        z <- z/max(z)
        z0 <- 0 #min(z) - 20
        z <- rbind(z0, cbind(z0, z, z0), z0)
        x <- c(min(x) - 1e-10, x, max(x) + 1e-10)
        y <- c(min(y) - 1e-10, y, max(y) + 1e-10)
        clr=colorRampPalette(c("honeydew","lawngreen"))
        nbcol=50
        iclr=clr(nbcol)
        nrz <- nrow(z)
        ncz <- ncol(z)
        zfacet <- z[-1, -1]+z[-1, -ncz]+z[-nrz, -1]+z[-nrz, -ncz]
        facetcol <- cut(zfacet, nbcol)
        fill <- matrix(iclr[facetcol],nr=nrow(z)-1,nc=ncol(z)-1)
        fill[ , i2 <- c(1,ncol(fill))] <- "white"
        fill[i1 <- c(1,nrow(fill)) , ] <- "white"
        par(bg = "transparent")
        persp(x, y, z, theta = 35, phi = 25, col = fill, expand=5, 
              shade=0.75,ltheta=45 , scale = FALSE, axes = TRUE, d=1,  
              xlab="Year",ylab="Age",zlab="Selectivity", 
              ticktype="detailed", ...)
      }
      ix <- 1:length(yr)
      plot.sel(yr, age, exp(log_sel[log_sel[,1]==gear,-1]),main=paste("Gear", gear))
    })
    if(gear==1){ # WARNING - assumes gears laid out a certain way !!!
      saveFig("fig.time.varying.comm.sel")
    }else if(gear==2){
      saveFig("fig.time.varying.surv.sel")
    }
    par(op)
  }else{
    cat(paste("Error plotting 'Time Varying Selectivity' figure.  There are no MPD outputs loaded for scenario",scenario,"\n\n"))
  }
}

plotRuntimeStats <- function(type=0,ylab=""){
  # plots runtime stats for all scenarios.
  # assumes all scenarios have the same number of projected years
  # assumes the opList structure is used.
  # types:
  # 1 = objFun, 2 = max gradient, 3 = number of function evals, 4 = hangcode, any other value = exit code
  try(dev.off(),silent=T)
  if(type==1 | type==2 | type==3){ # use PLOTBUBBLES from PBSModelling for these ones
    if(type==1){
      # Objective function values
      # GREEN means value is positive GOOD
      # RED means a bad objective function value, i.e. returned -1,#IND (which is represented as 0.0 from GrMPE)
      dat <- abs(opList[[1]][[4]]$ObjectiveFunction)
      for(scenario in 2:length(opList)){
        if(opList[[scenario]][[6]]){ # if MPD results are loaded
          dat <- rbind(dat,opList[[scenario]][[4]]$ObjectiveFunction)
        }else{
          dat <- rbind(dat,NA)
          cat(paste("Error plotting 'Runtime Statistics' figure.  There are no MPD outputs loaded for scenario",scenario,", it is NULL on the plot.\n\n"))
        }
      }
      dat <- t(dat)
      colnames(dat) <- 1:length(opList)
      rownames(dat) <- ""
      plotBubbles(dat,dnam=F,cpro=F,ylab=ylab,clrs=c("green","red","black"),xaxt='n',yaxt='n')
      text(1:length(opList),1.04,dat,srt=-45,adj=1)
      text(1:length(opList),1,1:length(opList))
      title("Objective function values")
    }else if(type==2){
      # Maximum Gradient values
      # GREEN represents a good gradient, i.e. one that is smaller than .maxGrad
      # RED represents anything greater than .MAXGRAD
      dat <- abs(opList[[1]][[4]]$MaxGrad)
      for(scenario in 2:length(opList)){
        if(opList[[scenario]][[6]]){ # if MPD results are loaded
          dat <- rbind(dat,opList[[scenario]][[4]]$MaxGrad)
        }else{
          dat <- rbind(dat,NA)
          cat(paste("Error plotting 'Runtime Statistics' figure.  There are no MPD outputs loaded for scenario",scenario,", it is NULL on the plot.\n\n"))
        }
      }
      dat <- t(dat)
      colnames(dat) <- 1:length(opList)
      rownames(dat) <- ""
      dat <- ifelse(dat>.MAXGRAD,0,dat)
      dat <- ifelse(dat<.MAXGRAD,dat,-dat)
      plotBubbles(dat,dnam=F,cpro=F,ylab=ylab,clrs=c("green","red","red"),xaxt='n',yaxt='n')
      text(1:length(opList),1.04,dat,srt=-45,adj=1)
      text(1:length(opList),1,1:length(opList))
      title(paste("Maximum gradient values (<",.MAXGRAD,")"))
    }else if(type==3){
      # Number of function evaluations
      # GREEN means the number of function evaluations was greater than .FUNEVALS
      # RED means the number of function evaluations was less than .FUNEVALS
      dat <- opList[[1]][[4]]$nf
      for(scenario in 2:length(opList)){
        if(opList[[scenario]][[6]]){ # if MPD results are loaded
          dat <- rbind(dat,opList[[scenario]][[4]]$nf)
        }else{
          dat <- rbind(dat,NA)
          cat(paste("Error plotting 'Runtime Statistics' figure.  There are no MPD outputs loaded for scenario",scenario,", it is NULL on the plot.\n\n"))
        }
      }
      dat <- t(dat)
      colnames(dat) <- 1:length(opList)
      rownames(dat) <- ""
      dat <- ifelse(dat<.FUNEVALS,-dat,dat)
      plotBubbles(dat,dnam=F,cpro=F,ylab=ylab,clrs=c("green","red","black"),xaxt='n',yaxt='n')
      text(1:length(opList),1.04,dat,srt=-45,adj=1)
      text(1:length(opList),1,1:length(opList))
      title(paste("Number of function evaluations (>",.FUNEVALS,")"))
    }
  }else{
  # NOT USING PLOTBUBBLES!!
    if(type==4){
      # Hang codes
      plotcharCol <- ifelse(opList[[1]][[4]]$HangCode==1,"red","green")
      # GREEN means no error condition
      # RED means no improvement in function value when 10th to last value compared with
      #     current value.
    }else{
      # Exit codes
      plotcharCol <- ifelse(opList[[1]][[4]]$ExitCode==1,"green","red")
      # GREEN for normal exit - i.e. all derivatives satisfy conditions
      # RED for problem with the initial estimate for the Hessian matrix.
      #     - The hessian matrix must be positive definite
      plotcharCol <- ifelse(opList[[1]][[4]]$ExitCode==2,"blue",plotcharCol)
      # BLUE for problem with the derivatives, either:
      # a) There is an error in the derivatives or
      # b) function does not decrease in direction of search, perhaps due to numerical
      #    round off error, or too stringent a convergence criterion
      plotcharCol <- ifelse(opList[[1]][[4]]$ExitCode==3,"purple",plotcharCol)
      # PURPLE for Maximum number of function calls exceeded
    }
    #  par( oma=c(2,2,4,1), mar=c(3,3,3,1), mfrow=c(1,1) )
    plot(1,1,
         pch=.PCHCODE,
         xlab="Scenario",
         ylab=ylab,
         col=plotcharCol,
         xlim=c(1,length(opList)),
         ylim=c(1,1))
    for(scenario in 2:length(opList)){
      if(type==4){
        plotcharCol <- ifelse(opList[[scenario]][[4]]$HangCode==1,"red","green")
        title("Hang code values")
      }else{
        plotcharCol <- ifelse(opList[[scenario]][[4]]$ExitCode==1,"green","red")
        plotcharCol <- ifelse(opList[[scenario]][[4]]$ExitCode==2,"blue",plotcharCol)
        plotcharCol <- ifelse(opList[[scenario]][[4]]$ExitCode==3,"purple",plotcharCol)
        title("Exit code values")
      }
      points(scenario,1,pch=.PCHCODE,col=plotcharCol)
    }
  }
}

######~~~~~~~~Unused in 2012 so far~~~~~~~~~~~###################################

fig.mcmc.diagnostics	<-	function(){
  #This function runs diagnostic plots for the posterior samples
  op <- par(no.readonly=T)
  if(length(unique(A$mc$abar))==1){
    # abar was not estimated, so don't include it
    if(length(unique(A$mc$gbar))==1){
      # gbar was not estimated, so don't include it
      np <- c(1:5,9)
    }else{
      np <- c(1:6,9)
    }
  }else{
    np <- c(1:7,9)
  }

	xyplot(post.samp[,np])
	acfplot(post.samp[,1:np])
  saveFig("fig.mcmc.diagnostics")
	par(op)
}

fig.fmsy.steepness <- function(){
  op <- par(no.readanly=T)
	par(mfrow=c(2,2), mai=c(0.75,0.751,0.15,0.15), omi=c(0.35,0.45,0.15,0.15), las=1)
	fmc=A$mc[,2]
	hmc=A$mc[,11]
	CRmc=4*hmc/(1-hmc)
	Mmc=A$mc[,3]

	hist(fmc,prob=T,breaks=30,xlab="Fmsy",ylab="Posterior density Fmsy",main="",col="wheat", cex.axis=1.1, cex.lab=1.1)
	lines(density(fmc,na.rm=T),col="red")
	mtext(paste("Mean =",round(mean(fmc),2), " Median =", round(median(fmc),2), "SD =", round(sd(fmc),2)), side=3, line=-0.5, cex=0.8)

	hist(hmc,prob=T,breaks=30,xlab="h",ylab="Posterior density Steepness",main="", xlim=c(0.2,1),col="wheat", cex.axis=1.1, cex.lab=1.1)
	lines(density(hmc,na.rm=T),col="red")
	mtext(paste("Mean =",round(mean(hmc),2), " Median =", round(median(hmc),2), "SD =", round(sd(hmc),2)), side=3, line=-0.5, cex=0.8)

	hist(CRmc,prob=T,breaks=30,xlab="CR",ylab="Posterior density Compensation Ratio",main="",col="wheat", cex.axis=1.1, cex.lab=1.1)
	lines(density(CRmc,na.rm=T),col="red")
	mtext(paste("Mean =",round(mean(CRmc),2), " Median =", round(median(CRmc),2), "SD =", round(sd(CRmc),2)), side=3, line=-0.5, cex=0.8)

	hist(Mmc,prob=T,breaks=30,xlab="M",ylab="Posterior density M",main="",col="wheat", cex.axis=1.1, cex.lab=1.1)
	lines(density(Mmc,na.rm=T),col="red")
	mtext(paste("Mean =",round(mean(Mmc),2), " Median =", round(median(Mmc),2), "SD =", round(sd(Mmc),2)), side=3, line=-0.5, cex=0.8)
  saveFig("fig.fmsy.steepness")

  par(op)
}

fig.spr.vs.management.target <- function(){
	#The relative spawning potential ratio (1-spr)/(1-spr.at.msy)
	op	<- par(no.readonly=T)
	spr <- A$mc.f40spr #read.table("tinss.f40spr",h=F)
	post.spr <- as.data.frame(window(mcmc(spr),start=Burn,thin=thin))
	sprci <- apply(post.spr,2,quantile,probs=c(0.025,0.5,0.975))
	matplot(A$yr,t(sprci),type="l",col=c(2,1,2),lty=c(3,1,3), lwd=2, pch=c(-1, 0, 1),ylim=c(0,max(sprci))
		,xlab="Year",ylab="SPR")
	abline(h=1)
	text(1980, 1, "Management target", pos=3)
	par(op)
}

fig.yields.4panel <- function(A,type){
	#plot the equilibrium yield curves
	op <- par(no.readonly=T)
	par(mfcol=c(2,2))
	#A$equil comes from the TINSS.rep file
	fe <- A$equil[, 1]
	ye <- A$equil[, 2]
	de <- A$equil[, 3]
	spr <- A$equil[, 4]

	plot(fe, ye, type="l", xlab="Fishing mortality (Fe)", ylab="Equilibrium yield"); gletter(1)
	plot(de, ye, type="l", xlab="Spawning depletion", ylab="Equilibrium yield", lty=2, col=2);gletter(2)
	plot(spr,ye, type="l", xlab="Spawning potential ratio", ylab="Equilibrium yield", lty=3, col=3);gletter(3)
	matplot(cbind(fe, de, spr), ye/max(ye)*100, type="l",xlab="Fe, depletion,  SPR",  ylab="Relative equilibrium yield")
	gletter(4)
  saveFig("fig.yields.4panel")
	par(op)
  
}

fig.yield.depletion.relrecuitment.spr <- function(){
	#Relationship between fishing mortlaity ~ yield,  recruitment,  SBe,  SPR
	op <- par(no.readonly=T)
	par(mfcol=c(2,2))
	fe <- A$equil[, 1]
	ye <- A$equil[, 2]
	de <- A$equil[, 3]
	spr <- A$equil[, 4]
	re <- A$equil[, 5]  
	ix <- c(min(which(ye==max(ye))),min(which(de<=0.4)) , min(which(spr<=0.4)))
	plot(fe, ye, type="l",xlab="", ylab="Equilibrium yield (million mt)", lwd=2)
	segments(fe[ix],0,fe[ix],ye[ix],lty=c(1, 2, 3))
	segments(0,ye[ix],fe[ix],ye[ix],lty=c(1, 2, 3))
	
	re <- re/re[1]
	plot(fe, re, type="l",xlab="", ylab="Relative recruitment", lwd=2) 
	segments(fe[ix],0,fe[ix],re[ix],lty=c(1, 2, 3)) 
	segments(0,re[ix],fe[ix],re[ix],lty=c(1, 2, 3)) 
	
	plot(fe, de, type="l",xlab="", ylab="Spawning depletion", lwd=2)
	segments(fe[ix],0,fe[ix],de[ix],lty=c(1, 2, 3))
	segments(0,de[ix],fe[ix],de[ix],lty=c(1, 2, 3))

	plot(fe, spr, type="l",xlab="", ylab="Spawning Potential Ratio", ylim=c(0, 1), lwd=2)
	segments(fe[ix],0,fe[ix],spr[ix],lty=c(1, 2, 3))
	segments(0,spr[ix],fe[ix],spr[ix],lty=c(1, 2, 3))   
	
	legend("topright", c("MSY", "SB40", "SPR40"), lty=1:3, bty="n")
	
	mtext("Equilibrium fishing mortality rate", 1, outer=T, line=-1)                               
  saveFig("fig.yield.depletion.relrecruitment.spr")
	par(op)
}

gweke.chain <- function (x, frac1 = 0.1, frac2 = 0.5, nbins = Nbin, pvalue = 0.05, auto.layout = TRUE, ...){
    x <- as.mcmc.list(x)
    oldpar <- NULL
    on.exit(par(oldpar))
    if (auto.layout) 
        oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), 
            Nparms = nvar(x)))
    ystart <- seq(from = start(x), to = (start(x) + end(x))/2, 
        length = nbins)
    if (is.R()) 
        gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 
            dimnames = c(ystart, varnames(x), chanames(x)))
    else gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)), 
        dimnames = list(ystart, varnames(x), chanames(x)))
    for (n in 1:length(ystart)) {
        geweke.out <- geweke.diag(window(x, start = ystart[n]), 
            frac1 = frac1, frac2 = frac2)
        for (k in 1:nchain(x)) gcd[n, , k] <- geweke.out[[k]]$z
    }
    climit <- qnorm(1 - pvalue/2)
    for (k in 1:nchain(x)) for (j in 1:nvar(x)) {
        ylimit <- max(c(climit, abs(gcd[, j, k])))
        plot(ystart, gcd[, j, k], type = "p", xlab = "First iteration in segment", 
            ylab = "Z-score", pch = 4, ylim = c(-ylimit, ylimit), 
            ...)
        abline(h = c(climit, -climit), lty = 2)
        if (nchain(x) > 1) {
            title(main = paste(varnames(x, allow.null = FALSE)[j], 
                " (", chanames(x, allow.null = FALSE)[k], ")", 
                sep = ""))
        }
        else {
            title(main = paste(varnames(x, allow.null = FALSE)[j], 
                sep = ""))
        }
        if (k == 1 && j == 1) 
            oldpar <- c(oldpar, par(ask = ask))
    }
 }
